An alternative to grids and glasses: 
Quaquaversal pre-initial conditions for N-body simulations 
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ABSTRACT 

N-body simulations sample their initial conditions on an initial particle distribution, which 
for cosmological simulations is usually a glass or grid, whilst a Poisson distribution is used for 
galaxy models, spherical collapse etc. These pre-initial conditions have inherent correlations, 
noise due to discreteness and preferential alignments, whilst the glass distribution is poorly 
defined and computationally expensive to construct. We present a novel particle distribution 
which can be useful as a pre-initial condition for N-body simulations, using a simple construction 
based on a "quaquaversal" tiling of space. This distribution has little preferred orientation (i.e. 
is statistically isotropic), has a rapidly vanishing large scale power-spectrum (P(fc) '-^ k^), and 
is trivial to create. It should be particularly useful for warm dark matter and cold collapse 
simulations. 



Subject headings: 

1. INTRODUCTION 

Numerical simulations have become a very pow- 
erful tool for investigating non-linear gravitational 
phenomenon, such as understanding the evolution 
and properties of cosmological structures. Since 
the simulations contain a rapidly increasing num- 
ber of particles, and probe structures on ever 
smaller scales, it is timely to address some of the 
fundamental aspects of the initial conditions used 
in these simulations. 

One such aspect is that the cosmological stan- 
dard model assumes that the early Universe is 
statistically isotropic. This is in contrast with 
the usual choice of pre-initial conditions for N- 
body simulations, most often given by placing 
particles on a uniform grid, which is intrinsi- 
cally anisotropic. While it has been shown that 
this anisotropy can produce non-physical effects in 



simulations, such effects are difficult to quantify. 
For cold dark matter simulations it is thought that 
the physically relevant correlations quickly grow 
and dominate over fluctuations due to discrete- 
ness. The same is not true for warm dark matter 
simulations for example. 

To address this and related questions numer- 
ically it is useful to have alternative pre-initial 
conditions (PrelCs). We construct a novel pre- 
initial condition by making use of a tiling of 
three dimensional space called the quaquaversal 
tiling (|Conwav fc RadinI 119981 ). We will show 
that this new particle distribution, which is sta- 
tistically isotropic (i.e. has little intrinsic direc- 
tionality), has mass fluctuations which decay as 
rapidly as in a grid. At the same time it is a 
deterministic structure, and can be trivially gen- 
erated. We present a C-code for the generation 
of these structures, which can be downloaded from 
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http : / /krone . phys ik . unizh . ch/ "hansen/ qua/ 



2. EXISTING PRE-INITIAL CONDI- 
TIONS 

Let us briefly recall the steps in any cosmologi- 
cal simulation. 1) First one chooses PrelCs, which 
most often is a regular grid, i.e. the particles are 
placed on a lattice. 2) One then imprints a power- 
spectrum onto these particles, by applying an ap- 
propriate displacement field specified through the 
Zeldovich approximation. 3) Then one runs the 
cosmological code, taking care of the many numer- 
ical issues related to convergence (softening, time 
stepping etc). The present paper focuses solely on 
the first step, namely the setting up of the pre- 
initial condition. 

There are three PrelCs which are regularly used 
in the literature. The first, which is the standard 
choice, places the particles on a regular grid. This 
is a very well tested method, which most likely 
produces the corr ect growth of long ran ge large 



scale fluctuations (jEfstathiou et al 
it 



19851) . How- 
ever, It IS unknown how much the grid affects 
small scale structures. One explicit example in 
which such effects have been observed to be im- 
portant are in simulations with warm dark mat- 
ter (WDM). In the WDM case the thermal mo- 
tion induces a free streaming, which e rases struc- 
tures on small scales. In a paper by iBode et al 



([20011) it was suggested that WDM might have the 
novel property of creating structures along a cos- 
mic web, b elow the cut - off fre quency of the power 
spectrum. iBode et al.l (|200lh took great care m 
testing for a large range of known numerical issues, 
and concluded that the effect observed was real 
and physical. It was later shown that these small 
scale structures were spurious . In the paper by 
Gotz fc Sommer-LarserJ (|2003l ) two almost iden- 
tical simulations were performed, differing only in 
the pre-initial conditions: one was a grid, the other 
was a glass (we will discuss glass initial conditions 
below). It was shown from the results of these 
simulations ( Gotz fc Sommer-LarsenI 20031) . that 
the conclusions reached in lBode et al.l (j200ll ) were 
incorrect precisely because of effects coming from 
the pre-initial particle distribution: the bead-like 
structures along the filaments observed were vir- 
tually absent in the simulation with a pre-initial 
glass. We emphasise that this is just an illus- 



tration of how difficult WDM simulations can be. 
Similar difficulties have been discussed in the con- 
text of cold dark matter (CDM) simulations e.g. 
start ing from specific config urations dMelott et al.l 
1997 ). and at early times ( Jovce et al.l 20051 ). it 



remains unclear how important such effects are in 
real simulations. 

The second PrelCs which are often used is the 
glass. The idea is to evolve a set of particles, ini- 
tially randomly distributed in a box, under nega- 
tive gravity, i.e. Coulomb forces, until one reaches 
a configuration in whic h the force o n each parti- 
cle is extremely small ([White I1994I 1. While this 
appears simple at first sight, there are a range 
of well known practical problems. Firstly, start- 
ing from the random configuration, the particles 
stream towards a lower potential, but gain kinetic 
energy which makes them oscillate about the min- 
imum of their local potential. To reduce the as- 
sociated Poisson noise, one needs to damp these 
velocities. If one does so by reducing the parti- 
cle velocities at a given time, then a large frac- 
tion of the particles will lie far from the mini- 
mum of the potential and little is gained. If one 
applies a more continuous damping of the parti- 
cle velocities, then either the small scale fluctua- 
tions are erased (if the damping is large), or the 
glass takes an unreasonably long time to create 
(if the damping is small). If one uses the method 
of simulated annealing, repeatedly heating up and 
cooling down the system, the creation of the glass 
becomes computationally very expensive. A ma- 
jor difficulty is that the final configuration is not 
unique, and indeed that one does not have a well 
defined criterion for determining when an optimal 
configuratio n has been reached. W e note in this 
respect fsee iGabrielli et al.l ([20031 )) that the sys- 



tem which is simulated is just a damped variant 
of the "one component plasma" (i.e. point par- 
ticles inter acting thro ugh Coulomb forces), which 
is known ([Carnll96lh to undergo a transition to 
a body centered cubic lattice configuration at low 
temperature. The glass is presumably a transient 
to such a grid type configuration, but it is not 
known what the relevant time scales are. 

A third possibility, and one that is frequently 
used to create equilibrium and non-equilibrium 
halo models is to simply select random positions 
for particles. In this case there is Poisson noise at 
all scales. If one tries to simulate a "cold collapse" 
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using such initial conditions then the growth of 
small scale structures can clearly be seen which 
will affect the virialisation of the final structure. 

In summary the grid PrelCs are trivial to cre- 
ate and has vanishing power-spectrum (below the 
Nyquist frequency). It has, however, a strong 
orientation and power on the scale of the grid 
spacing. The glass PrelCs have, in principle, lit- 
tle preferred orientation, and has a rapidly de- 
creasing power spe ctrum of density fl uctuations, 
with P{k) ~ (|Smith et all l2003l ). at large 
scales. These latter configurations are, however, 
not clearly defined and they are computationally 
expensive to create. Note, however, that the glass 
PrelCs still suffer from large scale anisotropy as 
long as periodic boundary conditions are used. 
Such large scale anisotropics are evidently also 
present in the new PreIC to be discussed below. 
The random PreIC has no orientation, but it has 
significant intrinsic power on all scales {P{k) = 
constant)). We will now present a novel PreIC 
which is clearly defined and easy to generate, has 
a large scale power spectrum vanishing as in the 
glass configurations {P{k) ~ fc^), and has no pre- 
ferred orientation. 




Fig. 1. — The quaquaversal tiling after one step. 
The thick lines (coloured) show the rotated smaller 
triangles, which are rotated by 27r/3 and tt/2 re- 
spectively. 



3. CONSTRUCTING A QUAQUAVER- 
SAL TILING 



space, based on a triangular prism that is re- 
peatedly rotated about orthogonal axes by angles 
27r/3 and 7r/2. 

The principle of our construction of the PrelCs 
is simple. The quaquaversal tiling defines a divi- 
sion of space into equal volume cells. The parti- 
cle distribution obtained by assigning a particle of 
the same mass to each cell is highly uniform: the 
only source of fluctuations is the redistribution of 
the mass at the scale of the cell. In deed a well 
known argumen t, due to Zcldovich (IZeldovich 



19651 : IZeldovich fc Novikov 1983). shows that such 



The quaquaversal tiling (jConwav fc Radin 



19981 ) is a hierarchical tiling of 3 dimensional 



a local mass conservation constraint should lead to 
a power spectrum with the behaviour P{k) ^ k^ 
at small k. If one adds the further constraint that 
the centre of mass is locally conserved one expects 
to obtain P{k) ~ fc^ at small k. Here we will 
impose this constraint by placing the particle in 
each tile at its centre of mass. Furthermore, the 
tiling has little preferred orientation, a desirable 
property which will be inherited by the particle 
distribution. 

Consider a triangular tile, made from a 1, ^3, 2 
right-angle triangle, with depth 1/2. This tile can 
be decomposed into 4 identical tiles, all with ex- 
actly the same properties as the original "parent" 
tile, by placing 3 lines from the center of the long 
side (length 2), to the centers of the other 2 sides 
(length 1 and \/3) and to the right angle. Then 
all lines are extended in depth (in the 3rd dimen- 
sion). One can now choose one of two possibilities, 
either to rotate the two triangles at the short axis 
by 27r/3 about an axis in the depth dimension, or 
to rotate the two triangles touching the right an- 
gle by 7r/2 about an axis in the y-dimension (see 
figure [T]). Finally, one places two tiles next to each 
other, each wit h their choice of rotation . For fur- 
ther detail, see lConwav fc RadinI ( 19981 ). 

After this first tiling, we are left with 8 trian- 
gles, each identical to the original parent trian- 
gle, but a factor 8 smaller in volume. (The angles 
of the individual triangles are the same, only the 
sides are smaller by a factor 2). We can now re- 
peat this process again for each triangle, giving 
us 64 identical triangles. After N such iterations 
we will have 8''^ identical triangles, with the same 
properties as the original parent triangle. 

To obtain our configuration, we then place a 
particle in the center of volume of each trian- 
gle. Finally we place two parent tiles on top 
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of each other, to form a rectangular box with 
sides 1,1, V3- This final distribution, which is our 
new PreIC to be used in simulations, thus con- 
tains 2 X 8^ particles. From now on, we will 
refer to this kind of particle distribution as a 
Q-set. For simulations where periodic boundary 
conditions are needed, one can only make Q-sets 
with 16, 1 28, ■ ■ ■ , 4.2M, 33.6M, ■ ■ ■ p articles. It is 
proved in Conwav fc RadinI ( 19981 ) that, in the 
limit of an infinite number of iterations, the orien- 
tations of the tiles are essentially random (uniform 
in S0(3)). We infer that our distribution (with fi- 
nite, but large, N) will have little directionality 
(see discussion in section 14. 4|) . 

Naturally one can imagine similar constructions 
based on other tilings. In this paper, however, we 
will focus on this simple quaquaversal structure. 

3.1. COMPUTER CODE 

We will briefly describe the idea behind this 
code. Each triangle is uniquely defined through 
the definition of the spatial position (in 3-d space) 
of 4 corners (3 corners would suffice, of course). 
Thus, the original parent tile is defined through 
4 different 3— vectors. This is expressed through 
one 12-vector, V. The 8 sub-tiles in the next level 
of tiling, n + 1, are defined through a constant 
matrix, M , applied to each vector at level n such 
that 



V,{n + l)=MV{n), 



(1) 



where i = 1, • • • , 8 define the 8 new sub-tiles. Thus 
the constant matrix Ai has 128 entries, and it 
uniquely defines each tiling step. 

The code recursively applies M to the vectors 
(it recursively tiles each subtile) , until the desired 
level, N, has been reached. Then a particle is 
placed in the center of volume of that triangle, 
and that point is written to a file. Constructed 
in this way the code is very fast and requires very 
little memory. 

In our implementation the matrix M is con- 
stant. This implies that the specific rotations by 
fractions of tt are identical at each level of refine- 
ment. It would be possible to generalize the pro- 
cedure, by allowing the rotations to differ at each 
level, e.g. sometimes rotate by 47r/3 instead of 
27r/3. In that way one could achieve a slightly 
higher level of isotropization. 



We will now consider the properties of the par- 
ticle distribution in this rectangular box, in terms 
of mass variance and power spectrum. 

4. STATISTICAL PROPERTIES 
4.1. Mass Variance a\.j{R) 

Let us first analyse the amplitude of mass fluc- 
tuations in a sphere of radius R with respect to 
the average mass. If M{R) is the mass (for a dis- 
crete distribution, the number of particles) inside 
a sphere of radius R, the normalised mass variance 
is defined as 



(i?) = 



{MjHf) - {M{R)f 
{M{R)y 



(2) 



where the brackets indicate an ensemble average. 
For a distribution like ours (or e.g. a grid) one can 
define the ensemble average as the average over 
random positions of the initial box. Such an en- 
semble average definition is equivalent to a spatial 
average, with the mass variance then given as the 
infinite volume limit of the estimator defined be- 
low. 

For our mass variance calculations we have used 
the simple estimator 



•'j\/,cst 



(i?) 



1 



{Nr] 



E 



N,-l 



,(3) 



where N.i{R) is the number of particles in the ith 
of Ns randomly thrown spheres, constrained to be 
inside the sample volume. (Nr) is the mean num- 
ber of particles in such a sphere, given exactly by 
(Nr) — where v is the volume per particle 

(i.e. the volume of a single tile in the Q-set). 

We apply this estimator to a grid, a Pois- 
son distribution, a typical ACDM initial condition 
(z ~ 70), and a Q-set. We have used 128'^ parti- 
cles in a cube of side unity for the first three, while 
the Q-set employed seven tiling levels, which gives 
2 X 8"^ = 2 X 128^ particles in a rectangular box. 
The dimension of the Q-set box is chosen so that 
the mean density is equal to that of the other dis- 
tributions. This is a convenient choice for compar- 
ison of the results, as it is the mean particle density 
riQ which fixes the asymptotic level of the Poisson 
variance at small scales in any point distribution. 



with <7jj{R ^ 0) 



where Vs = 4ttR^/3 
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It is defined as 
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P(k)= lim -(|<5(k)|^ 

V^oo V 



(4) 



Fig. 2. — The mass variance as function of ra- 
dius. The correct (analytical) behaviour of the test 
structures are found. The Poisson (red, dotted) 
has ^ R^^, the grid (black, long-dashed) has 
cr^ ~ turning into for distances much be- 
low the interparticle distance. The ACDM (green, 
dashed) has the correct behaviour until scales 
where the interparticle distance approaches the 
grid size. The Q-set (blue, solid line) has the same 
behaviour as a grid, namely cr^ ~ on all scales 
down to the Poissonian turn over. 



(see e.g. Gabrielli. Jovce fc Svlos Labinil ( 2002f )). 
Several tests confirm that our estimator calculates 
the correct spatial properties (see further discus- 
sion below) . We see (figure [2]) that the mass vari- 
ance of a Q-set has the same behaviour as that 
of a grid with ^ above the interpar- 

ticle distance. The mean interparticle distance, 
A = 1/128, for the structures with 128^ particles 
is also shown in the figure. Note that this is ex- 
actly the length of the shortest side of a tile in the 
level seven quaquaversal tiling used for the Q-set 
considered here. 

4.2. Power Spectrum P{k) 

The power spectrum is the primary statistical 
tool used to characterise fluctuations in cosmology. 



where 6(h) is the Fourier transform of the density 
fluctuation fleld (5(x) = (p(x) — /9o)/po- In the case 
of a discrete distribution (i.e. of point particles) 
these quantities simply become 

V v-^ 

■ 2^ (5d(x - Xp) - 1 , (5-a) 



<5(x) = 
^(k) = 



iVp 



-ik-Xv, 



(k ^ 0) , (5-b) 



where Xp is the location of each particle and (5d is 
the Dirac delta function. 

To estimate the power spectrum we have used 
the "brute force" method i.e. we calculate it di- 
rectly from the formula one obtains by substi- 
tuting Eq. (j5-bp in Eq. Q, without the infinite 
volume limit and ensemble average. Our finite 
volume V is thus a rectangular box with sides 
Li, and we assume periodic boundary conditions 
so that k in the Fourier sums take the values 
k = 2tt (rix/ Lx,ny/ Ly^Hz/ Lz) where are in- 
tegers. We obtain P{k) = P(|k|) by ave raging 
over a bin of finite width around k — |k| ( Sirkd 
20051 ). It is important to take care in the inter- 



pretation of the large scale (i.e. small fc) modes, 
which will be affected both by under-sampling and 
contaminated by the boundary conditions (which 
systematically suppress power at small fc). With 
this simple estimator, however, we do not have 
to worry about the effect of assignment function, 
which typi cally causes problems in EFT methods 

As for the mass variance we calculate the power 
spectrum for a grid, a Poisson distribution, a typ- 
ical ACDM initial condition, and a Q-set. This 
time we use a smaller number of particles in order 
to facilitate the more computationally demanding 
procedure of calculating P(fc): they all have 32'^ 
particles except for the Q-set which has N = 5, 
and therefore again has twice as many points 
(32'^ — 8^). Likewise for the mass variance, this 
makes the the Poisson level {V/Np) the same for 
all the distributions considered. 

We see in figure [3] that the Q-set has the antic- 
ipated P{k) ~ fc^ behaviour up to the wave vector 
for the mean interparticle separation, fcA — 27r/A, 
where it flattens on average, to the Poisson level. 
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We have used a very fine binning in order to cap- 
ture this characteristic slope as well as the peaks 
arising from typical interparticle distances. We 
note that peaks appear at frequencies of 2^"fc, 
where n has integer values. It can be shown (see 
Radi that in hierarchical structures, a 

feature in real space between two radii, ri and r2, 
also must appear between the radii kti and Kr2, 
where k = 2 for a quaquaversal structure. 

Consistent with this interpretation wc note that 
adding random perturbations to the particles, 
hence creating a shuffled Q-set, gives P{k) ~ 
with peaks of diminished amplitude. 

A comparison between the different distribu- 
tions can be seen in figure ID Note that the grid 
has zero power up where it is strongly peaked 
and that all distributions, on average, reach the 
Poisson level at this wave number. 






Fig. 3. — Detailed plot of the power spectrum for 
a the Q-set. We have used a very fine binning to 
capture the slope as well as the characteristic 
peaks of this PrelC. 



4.3. Relations betvifeen <t'Ij{R) and P{k) 

a^(i?) and P{k) are related by the standard ex- 
pression 



1 



+ OC' 



P{k)W^{kr)k^dk, (6) 



Fig. 4. — The power spectra for a Q-set (solid, 
blue), a grid (short-dashed, black), a Poisson dis- 
tribution (dotted, red) and a ACDM initial con- 
dition (long-dashed, green) . The Q-set shows a 
fc^ behaviour up to the Poisson level. The power 
spectrums has characteristic peaks appearing, as 
discussed in the text, at fc/2, A;/4, fc/8 . . . due to 
the hierarchical properties of the Q-sct. 

where W'^{kr) is the Fourier transform of the 
spherical top hat window function. 

By studying Eq. ([6]) for a po wer spectrum of the 

form P(fc 0) ^ A;" one finds (jGabrielli. Joyce fc Svlos Labini 
2OO2I ) that for i? ^ 00 



'M 



(R) 



if n < 1 
log(i?)/i?4 if n = 1 
1/R* if n > 1 



(7) 



Two particular and simple examples which are 
useful reference points against which to gauge new 
distributions are: 



• Poisson: a 
const 



(R) 



R 



-3 



Shuflled Lattice: afjiR) 



P{k) ^ V/Np 



P{k) - P 



The numerical results for the mass variance 
and power spectrum in the previous two sec- 
tions are all in line with these analytic re- 
sults, notably our new Q-set has (j\j{R) ~ 
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and P{k) ^ fc^. This (large scale) behaviour 
makes the Q-set a member of a group of sys- 
tems which have been termed super- homogen ous 
( Gabrielh. Jovce fc Svlos Labinill2002l:lGabrieni etal 



2003 ) or hyperuniform (|Torauato fc Stillinger 
20031 ). Such distributions, defined by the prop- 
erty P{k ^ 0) = 0, are characterised in real 
space (cf. Eq. ^ above) by the asymptotic be- 
haviour of their variance, (J^{R) ^ with 
d < m < d -|- 1 where d is the spatial dimension. 
This quantity thus decays faster than in a Poisson 
point process (m = d), and in the cases we have 
considered attains the behaviour (m = d + 1). 
This is in fact the fastest possible decay of this 
quant ity for either point (.B eck 1 9871) or c ontin- 
uous (iGabrielh. Jovce k Svlos Labi3l2002l) mass 
distributions. We note that the glass PreIC also 
belong to this class, as it shares this same be- 
haviour of the variance and power spectrum as 
the Q-set we have introduced and analysed. 




Fig. 5. — 9 and (j) are the standard angles in spheri- 
cal coordinates, giving the orientations of a chosen 
reference vector in each tile of a level 7 Q-set (blue 
triangles), and the analagous quantities for a grid 
(red open circles). 

4.4. Isotropy 

One of the nice features of the Q-set is, as we 
have underlined, that it is isotropic. This result 



applies, however, in the limit of an infinite number 
of iterations of the algorithm we have described. It 
it interesting to quantify the degree of isotropy of 
a finite level tihng. To do so we show in Figure [5] a 
plot giving (blue triangles) the distinct directions 
defined by the tiles of the quaquaversal tiling used 
to construct the Q-set with N — 7. The directions 
are those of the shortest axis of each tile with re- 
spect to the orientation of the mother tile. Also 
shown (red open circles) is an analagous character- 
isation of a grid, which gives just the six orienta- 
tions of the vectors pointing to nearest neighbour 
sites. 

Figure [6] shows, on the other hand, the depen- 
dence of the number of distinct orientations, in the 
{6, (/))-plane, on the level of the tiling. The degree 
of isotropy grows very rapidly with increasing TV. 




Fig. 6. — Number of distinct orientations as a 
function of level of tiling. 



5. CONCLUSIONS 

We have studied an alternative to the stan- 
dard pre-initial conditions for N-body simulations 
of cosmological structures. The standard grid has 
strong orientations, and the glass is poorly defined 
and computationally expensive to create. We have 
therefore considered a particle distribution cre- 
ated starting from an equal volume tiling of space. 
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called the quaquaversal tiling. The particle distri- 
bution is trivial to create, has virtually no orien- 
tation (is statistically isotropic), and has rapidly 
vanishing large scale power-spectrum. We provide 
a C-code for the generation of these structures on 
[http : / /krone . phys ik . uniz h . ch/~hanse n/qua/[ 
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